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Abstract 

Uncertainties in L-band (1.4 GHz) radiative transfer modeling (RTM) 
affect the simulation of brightness temperatures (Tb) over land and the in- 
version of satellite-observed Tb into soil moisture retrievals. In particular, 
accurate estimates of the microwave soil roughness, vegetation opacity and 
scattering albedo for large-scale applications are difficult to obtain from held 
studies and often lack an uncertainty estimate. Here, a Markov Chain Monte 
Carlo (MCMC) simulation method is used to determine satellite-scale esti- 
mates of RTM parameters and their posterior uncertainty by minimizing 
the misfit between long-term averages and standard deviations of simulated 
and observed Tb at a range of incidence angles, at horizontal and vertical 
polarization, and for morning and evening overpasses. Tb simulations are 
generated with the Goddard Earth Observing System (GEOS-5) and con- 
fronted with Tb observations from the Soil Moisture Ocean Salinity (SMOS) 
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mission. The MCMC algorithm suggests that the relative uncertainty of 
the RTM parameter estimates is typically less than 25% of the maximum 
a posteriori density (MAP) parameter value. Furthermore, the actual root- 
mean-square-differences in long-term Tb averages and standard deviations 
are found consistent with the respective estimated total simulation and obser- 
vation error standard deviations of cr m = 3.1 K and a s = 2.4 K. It is also shown 
that the MAP parameter values estimated through MCMC simulation are 
in close agreement with those obtained with Particle Swarm Optimization 
(PSO). 

Keywords: 

radiative transfer modeling, brightness temperature, Bayesian parameter 
estimation, uncertainty, Markov Chain Monte Carlo simulation, SMOS 
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1. Introduction 


Uncertainties in radiative transfer modeling (RTM) affect the simula- 
tion of brightness temperatures (Tb) over land and the inversion of satellite- 
observed Tb to soil moisture retrievals. Quantification of these uncertainties 
is crucial to producing, validating and using passive microwave data, such 
as those obtained from the Soil Moisture Ocean Salinity (SMOS, Kerr et al. 
(2010)) and future Soil Moisture Active Passive (SMAP, Entekhabi et al. 
(2010)) missions. Yet, it is not particularly clear which RTM formulation 
and parameter values to use for large-scale applications. 

In the context of forward Tb simulation, several studies have analyzed 
the effect of different RTM formulations for the microwave roughness length, 
vegetation parameterization and soil dielectric model (Drusch et al., 2009; 
de Rosnay et al., 2009). The impact of parameter values and dynamic land 
surface variables as input to large-scale forward Tb simulations was demon- 
strated by, e.g., De Lannoy et al. (2013) and Balsamo et al. (2006), re- 
spectively. Similarly, soil moisture retrievals based on Tb observations are 
affected by the RTM formulation and parameter values (Crosson et al., 2005; 
Panciera et al., 2009; Konings et al., 2011; Parinussa et al., 2011), as well 
as by the choice of background and auxiliary fields, such as soil temperature 
and vegetation characteristics (Kerr et al., 2012; O’Neill et al., 2012). Col- 
lectively, these studies suggest that RTMs exhibit significant uncertainty and 
that the precise magnitude and impact of this uncertainty on large-scale Tb 
simulations and soil moisture retrievals remain unclear. 

Estimating the uncertainty of microwave RTM parameters is a major 
challenge, especially at larger spatial scales. Field experiments have pro- 
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vided RTM parameters values (de Rosnay et al., 2006; Grant et al., 2007; 
Panciera et al., 2009; Sabater et al., 2011), but mostly without uncertainty 
estimates. De Lannoy et al. (2013) derived global-scale RTM parameter val- 
ues and ad hoc uncertainty estimates using SMOS observations and Particle 
Swarm Optimization (PSO, Kennedy and Eberhart (1995)). PSO is espe- 
cially designed to find the optimal parameter values within a limited budget 
of function evaluations, but without recourse to estimating their underlying 
uncertainty. 

In this paper, we introduce a (Bayesian) Markov chain Monte Carlo 
(MCMC) simulation method to estimate the posterior RTM parameter dis- 
tribution. The DiffeRential Evolution Adaptive Metropolis (DREAM) algo- 
rithm is used with parallel direction and snooker sampling from past states 
(Vrugt et al., 2008, 2009; Laloy and Vrugt, 2012), referred to as DREAM( ZS ). 
Bayesian approaches such as DREAM( Z s) have many advantages over op- 
timization methods such as PSO. The explicit treatment and analysis of 
uncertainty help to understand which parts of the RTM model are well re- 
solved and which elements require further attention. Furthermore, a formal 
analysis of the residuals can be used to check the validity of our assump- 
tions about the residual error distributions and to discern whether reliable 
parameter values have been derived. 

The added value of obtaining posterior parameter distributions with Bayesian 
approaches, however, comes at an increased computational cost. Adequately 
sampling the posterior parameter distributions is too costly for global-scale 
operational applications that rely on evolving modeling systems in need of 
frequent re-calibrations, but can provide a valuable benchmark to verify re- 
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suits from simple parameter optimization algorithms, such as for example 
PSO. 

The goals of the present paper are thus to infer RTM parameters and 
their posterior uncertainty using a Bayesian method, and to study the as- 
sociated simulated Tb uncertainty. We are using the Goddard Earth Ob- 
serving System (GEOS-5) modeling framework that will be used to gener- 
ate the planned global SMAP Level 4 Surface and Root Zone Soil Moisture 
(L4_SM) data product through assimilation of SMAP Tb observations (Re- 
ichle et ah, 2012). As in De Lannoy et al. (2013), we focus on optimizing 
time-invariant RTM-parameters by minimizing climatological differences be- 
tween multi-angular, horizontally and vertically polarized Tb for morning 
and evening overpasses from SMOS observations and GEOS-5 simulations. 
The time-invariant optimized parameters will later be used in a data assimi- 
lation system (outside the scope of this paper), where state variables such as 
soil moisture and soil temperature will be updated in response to short-term 
variations in the observed Tb. 

To summarize, in this paper we apply MCMC simulation using multi- 
angular SMOS Tb observations to (i) verify if the maximum a posteriori 
density (MAP) parameter values derived from a converged posterior distri- 
bution with DREAM(zs) can be approximated using PSO, (ii) obtain reliable 
parameter uncertainty estimates, and (iii) quantify the magnitude of param- 
eter and other error sources in Tb simulations. The remainder of this paper 
is organized as follows. Section 2 summarizes the modeling system and the 
SMOS observations used in the present study. This is followed in section 3 
by a description of the DREAM( Z s) MCMC simulation method and PSO. 
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Section 3 also discusses several quantitative diagnostic metrics to analyze 
the simulated Tb uncertainty. Finally, this paper concludes in sections 4 and 

5 with a discussion of the results and conclusions. 

2. Observations and Model 

2.1. SMOS Tb Data 

Since its launch in November 2009, the SMOS mission provides global Tb 
data at a nominal spatial resolution of 43 km and with an equator overpass 
every 3 days. Here we use the multi-angular, full polarization Tb data from 
the period 1 July 2010 to 1 July 2012. Specifically, the data are extracted 
from the MIFLSCLF1C product, with processor version 504 for the years 2010 
and 2011, and version 551 from January 2012 onwards. Our previous study 
presented in De Lannoy et al. (2013) discusses in detail the various steps 
involved in the processing of the SMOS data. Most importantly, the data 
are screened extensively using both product-based data quality information 
and model-based quality control rules. Furthermore, the data are spatially 
mapped onto a 36 km Equal- Area Scalable Earth Grid (EASE) and binned 
per incidence angle. Consistent with our previous study, only a subset of 

6 incidence angles is used: 0=[32.5°, 37.5°, 42.5°, 47.5°, 52.5°, and 57.5°], 
where, for example, 32.5° represents the average of all data with incidence 
angles between 32° and 33°. 

For the purpose of estimating the microwave RTM parameters, long-term 
averages ( m 0 ) and standard deviations (s G ) of the SMOS data are computed 
separately for each of the 6 incidence angles, 2 polarizations (horizontal H 
and vertical V), and 2 overpass times (ascending at 06:00h local time (LT), 
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descending at 18:00h LT). This results in a total of 48 “observations” per grid 
cell: 24 for the long-term average Tb and 24 for the long-term Tb standard 
deviation. Section 3 provides more details. 

2.2. GEOS- 5 Tb Modeling 

The modeling combines (i) land surface modeling with the Catchment 
land surface model (CLSM) and (ii) radiative transfer modeling with a tau- 
omega model to simulate long-term Tb averages and standard deviations. As 
in De Lannoy et al. (2013), the GEOS-5 CLSM (Koster et ah, 2000) is set up 
on the 36 km EASE grid and spun up prior to the SMOS observation period. 
Surface meteorological forcing data at a 1/2° x 2/3° spatial and hourly tem- 
poral resolution are taken from the Modern-Era Retrospective analysis for 
Research and Applications (MERRA, Rienecker et al. (2011)). The MERRA- 
precipitation is corrected with gauge-based precipitation from the National 
Oceanic and Atmospheric Administration (NOAA) Climate Prediction Cen- 
ter “Unified” (CPCU) product (Reichle, 2012). The model version is the 
same as that used for the MERRA-Land data product (Reichle et al., 2011), 
except for two changes that more closely align the model with the version 
that will ultimately be used for the SMAP L4_SM data product: (i) the sur- 
face soil moisture pertains to the top 5 cm surface layer (as opposed to the 
top 2 cm layer in MERRA-Land), and (ii) a preliminary version of updated 
soil parameters from a forthcoming version of GEOS-5 is used. 

The vegetation parameterization in CLSM uses 8 default vegetation classes. 

For the RTM simulations, these classes are further refined into the 16 classes 
defined by the Moderate Resolution Imaging Spectroradiometer (500 m MOD12Q1 
V004) International Geosphere-Biosphere Programme (IGBP) land cover 


7 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 


classification (Loveland and Belward, 1997). Figure 1 shows the North Amer- 
ican study domain which covers 9 of the 16 IGBP vegetation classes. 

The soil moisture, soil temperature, vegetation water content, air temper- 
ature and climatological vegetation dynamics simulated with the prognostic 
CLSM are used as input to the diagnostic zero-order (tau-omega) microwave 
RTM to simulate L-band Tb. A short description of the RTM is given in Ap- 
pendix. In essence, the Tb is determined by the surface soil temperature and 
attenuated by dynamic and static soil and vegetation characteristics. The 
key model parameters that impact the rough surface reflectivity h (Eq. A. 3, 
Eq. A. 4), the scattering albedo a;, and vegetation optical depth r (Eq. A. 6) 
will be estimated using the multi-angular SMOS observations (section 3), 
where h is a function of soil moisture and r depends on the leaf area index 
(LAI). 

3. Methods 

3.1. Overview 

Keeping up with our previous work (De Lannoy et ah, 2013), the objec- 
tive of the parameter estimation is to minimize the difference between long- 
term (climatological) averages and standard deviations for multiple types of 
SMOS-observed and GEOS-5-modeled Tb. We purposely do not minimize 
differences in the time domain as the goal of the present paper is to derive 
parameter estimates that result in the smallest possible bias in the long-term 
simulation of Tb. Short-term differences between Tb observations and simu- 
lations will be exploited in future studies using sequential data assimilation. 
We estimate a time-invariant multi-dimensional parameter set (hereafter re- 
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ferred to as a) that determines climatological features of the simulated Tb. 
The parameters are optimized locally, i.e. , for each grid cell independently, 
and only for non-frozen land surface conditions as determined by the GEOS-5 
modeling system. 

Table 1 gives an overview of the parameters estimated in different exper- 
iments. All scenarios estimate the 5 most relevant RTM-parameters: h min , 
Ah = h max — h min , &#, A b = by — bu and u (according to the best scenario 
identified in De Lannoy et al. (2013)). Based on these time-invariant pa- 
rameters, time-variant values of /i, th and tv are computed, using dynamic 
information about soil moisture for h (Eq. A. 4) and LAI for r (Eq. A. 6). 
Time-averaged results for < h > and < r > are then presented, where < • > 
denotes the long-term time average. These RTM-parameters are estimated 
with either DREAM(zs) or PSO, hereafter referred to as scenarios D and P, 
respectively. The DREAM( Z s) analysis is further expanded to also include 
the residual Tb error statistics cr m and a s (scenario D^, discussed below). 
For each grid cell, we thus estimate 5 parameters for scenarios P and D, and 
7 for D a . 

To derive these parameters, we minimize per grid cell the climatological , 
or long-term, differences between 48 Tb observations and simulations. The 
2 x 24 observations consist of long-term Tb averages and Tb standard de- 
viations for the 24 combinations of 2 polarizations, 2 overpass times, and 6 
incidence angles. The errors in these observations are assumed to be indepen- 
dent, that is, we neglect correlations in instrument errors and errors between 
H- and V-polarized observations at identical incidence angles. Similarly, the 
simulation errors are assumed to be independent, even though correlation is 
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to be expected. Note that temporal correlations in the errors are of little 
concern because the observations are long-term averages and standard de- 
viations, and not measurements in the time domain (Wohling and Vrugt, 
2011). 

In keeping up with De Lannoy et al. (2013), the two years of historical 
SMOS data are divided into a calibration period (1 July 2011 - 1 July 2012) 
and an evaluation period (1 July 2010 - 1 July 2011). To ensure a meaningful 
calibration at each grid cell, we impose a minimum of 20 valid data points 
(Ni) per year to compute the long-term Tb average and standard deviation 
for a particular combination (z = 1, . . . , 24) of incidence angle, polarization 
and overpass time. The requirement of Ni > 20 is used for the calculation of 
evaluation statistics as well. 


3.2. Markov Chain Monte Carlo (MCMC) Sampling 

The Bayesian framework allows deriving posterior probabilities of param- 
eter estimates and model simulations, conditioned on errors in observations 
and simulations. The posterior probability distribution is computed by com- 
bining the observation likelihood p(m 01 s 0 \cx ) with a prior distribution p(cx): 


p(a|m 0 , s 0 ) = 


p(m 0 , s 0 \cx)p(cx) 


(i) 


J a p(m 0 , s 0 \ct)dct 

The observations consist of long-term averages (m ij0 E m 0 ) and standard 
deviations (s^ 0 E s Q ) of Tb for 24 different combinations of incidence angles, 
polarizations and overpass times (z = 1, . . . , 24). The denominator is a nor- 
malization factor and thus it suffices to maximize p(m 0 , s 0 \cx.)p(cx) to find 
the posterior distribution of a. In practice, it is difficult to solve this prob- 
lem analytically and we therefore resort to MCMC simulation to generate a 
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sample of the posterior target distribution. 

In this paper, the differential evolution adaptive metropolis (DREAM( Z s), 
Vrugt et al. (2008); Laloy and Vrugt (2012)) algorithm with sampling from 
past states is used to efficiently explore the posterior parameter distribution. 
This algorithm adaptively updates the scale and orientation of the proposal 
distribution during sampling, and is specifically designed to rapidly explore 
multi-dimensional target distributions. In DREAM( Z s), multiple chains are 
running in parallel and the update of a chain is determined from an external 
sample of points that collectively summarizes the search history of all the 
individual chains. The log-likelihood of the current and proposed parameter 
values are compared using the Metropolis selection rule. If the proposal is 
accepted, the chain moves to this new point, otherwise the chain remains 
at its current position. Diminishing adaptation of the external archive of 
samples ensures convergence to the exact posterior distribution. 

We assume a Gaussian prior for each of the individual parameters olq^ E 
ctQ. The mean and standard deviation of this multi-normal distribution p{cx) 
are derived from literature values that yield reasonable Tb simulations com- 
pared to SMOS Tb and are summarized in Table 1. Note that these values 
were referenced as ‘Lit2’ in De Lannoy et al. (2013). The prior mean for each 
individual parameter is given by a vegetation-dependent value a and the 
standard deviation is defined by a 2 aQ k = ( a max ,k - a m in,k) 2 / 12, using 

upper and lower bounds [a maXtk , a mimk ]. 

The following log-likelihood function is used to minimize the differences in 
long-term Tb averages and standard deviations between observations (m^ G , s ijG ) 
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and simulations (m^(a), Si(cx)): 



L s ,o (2) 


This formulation thus explicitly takes into consideration long-term biases in 
the Tb average (L mj0 [-]) and the Tb variability (L 5?0 [-]) and is derived from 
a classical Gaussian likelihood function: 


where and denote the (ensemble) standard deviations of the residual 


averages and standard deviations, respectively. 

3.3. Particle Swarm Optimization (PSO) 

The PSO algorithm (Kennedy and Eberhart, 1995) is a global search 
method that uses a dynamic swarm of particles to explore the parameter 
space. The best position of each individual particle (cognitive aspect) and 
of the entire swarm (social aspect) are used to guide the particles towards 
the optimal solution. The iterative swarm search is performed in several 
independent repetitions to account for sampling variability. 

The fitness of each parameter combination in the swarm is measured by 
an integrated Tost’ or ‘objective function’ J [-] that measures the distances 



(3) 


differences between the observed and simulated values of the long-term Tb 
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between the observed and simulated long-term Tb averages (J m>0 [-]) and 
standard deviations (J S:0 [-]). To make sure that the estimated parameter 
values honor the prior information (as used in DREAM( Z s)), we also include 


pected values (J a [-]). This results in the following definition of the objective 
function to be minimized: 


where N a signifies the number of simultaneously estimated parameters. This 
formulation is essentially similar to the definition of the posterior density 
used in DREAM( Z s). The main difference is that PSO handles the prior 
information of the parameters explicitly as penalty term J a in the objective 
function, whereas in DREAM( Z s), the prior parameter distribution is handled 
independently from the likelihood function by application of Bayes law. Both 
methods should thus find the same “best” parameter values. 

3.4- Likelihood , Objective Function and Algorithm Settings 

The design of the likelihood (L) or objective ( J) function for DREAM( Z s) 
and PSO warrants further discussion. As discussed above, we sample the cli- 
matological, or long-term, Tb averages and standard deviations over multiple 
incidence angles, polarizations and overpass times (that is, 2 x 24 observa- 
tions, i = 1,...,24) per location, rather than one observation at multiple 


a penalty term that quantifies deviations of the parameters from their ex- 
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time steps. The long-term Tb averages and standard deviations could also 
be interpreted as ‘summary statistics’ or ‘signatures’ of the system, and hence 
our approach has many elements in common with the diagnostic model eval- 
uation procedure presented in Vrugt and Sadegh (2013). 

The variables ay m and a^ s in Eq. 2 and Eq. 4 measure the (ensemble) 
standard deviation of the residual differences between the observed and sim- 
ulated long-term Tb averages and standard deviations, respectively, for each 
observation i. The residual errors are assumed to have a zero mean and in- 
clude both SMOS observation and simulation errors, due to e.g. inaccurate 
soil moisture, temperature or vegetation characteristics. These and 
statistics trade-off errors in the long-term Tb averages against those of the 
long-term Tb standard deviations (as well as deviations from the prior pa- 
rameter constraints). Since only one sample is available for each observation, 
it is impossible to estimate individual oy m - and cq ?s -values. Therefore, we de- 
fine and as a combination of a homoscedastic term (cr m , <j s ) and a 
tuning factor Wi to account for the robustness of the diagnosed long-term Tb 
averages and standard deviations, i.e. cr 2 m = and af s = w *cr 2 . The 

homoscedastic term is identical for all 24 observations and set to a default 
value of 1 K (De Lannoy et ah, 2013), or alternatively we estimate a m and a s 
jointly with the RTM parameters (see section 3). The weights are given by 
Wi = ;, where Ni denotes the number of data points in time that contribute 

to a particular long-term Tb average (or standard deviation), and N signifies 
the average number of time steps across all observations. These weights are 
typically close to 1 and assign somewhat more (less) weight to climatologi- 
cal differences that are based on more (fewer) individual data points in the 
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original 1-year data time series. 

A maximum of 12,000 log-likelihood function evaluations are performed 
with DREAM( Z s) using standard settings of the algorithmic variables. For 
PSO, we use the same algorithmic settings as De Lannoy et al. (2013), except 
a swarm size of 10 particles is used with a minimum of 10 and maximum of 100 
iterations. The search is terminated if the reduction of the objective function 
is smaller than IE-5 over the last 10 iterations. A total of 12 repetitions are 
performed, which results in a maximum of 12,000 function evaluations. 

3.5. Posterior Parameter Distribution 

The ‘optimal’ parameter values are defined as those with the maximal a 
posteriori density (MAP), i.e. with the largest value for L (Eq. 2, DREAM( Z s)) 
or smallest value for J (Eq. 4, PSO). Note that these MAP values are not nec- 
essarily identical to the posterior ensemble mean of the distribution derived 
with of DREAM(zs)- For the DREAM(zs) experiments, the last 25% of the 
MCMC chains (3,000 samples) are used to summarize parameter uncertainty 
by calculating the standard deviation of each individual parameter. To illus- 
trate this in more detail for one grid cell, consider Fig. 2a, which depicts the 
marginal distributions of the RTM parameters. We define the uncertainty as 
the ensemble standard deviation stdv[a] = a — a centralized around the en- 
semble mean d 7 , not around the MAP parameter value a map- The notation 
“ refers to the ensemble mean. Note that the standard deviation around 
the MAP estimate stdvMAp[] can be found as a function of the centralized 
standard deviation stdv[.\, i.e. stdvMAp[ -] 2 = stdv[.] 2 + A a (A a — stdv[.]), 
where A a = a — olmap is the difference between the ensemble mean and 
MAP parameter estimate. We found that, across the different experiments, 


15 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 


A a is either small or A a and stdv[.\ are of similar magnitude (not shown 
herein), so that stdvMAp[] ~ stdv[.]. 

3.6. Convergence 

‘Convergence’ can reflect accuracy (closeness to the actual optimum solu- 
tion) or precision (reduction of the prior uncertainty). The following hypothe- 
ses will be verified to assess the convergence of the DREAM( Z s) algorithm: 
(i) the Tb performance (accuracy) with posterior parameter estimates should 
be better than with prior parameter guesses (section 3.7), (ii) the posterior 
parameter uncertainty (section 3.5) and the corresponding uncertainty in Tb 
simulations (section 3.7) should be reduced compared to their counterparts 
derived from the prior parameter uncertainty, and (iii) the potential scale re- 
duction factor VR by Gelman and Rubin (1992) should be near 1 to inspire 
confidence that the different MCMC chains have converged to the appro- 
priate limiting distribution. The latter metric measures by which scale the 
posterior distribution will shrink as the number of MCMC iterations would 
go to infinity. 

3. 1. Tb Performance and Ensemble Verification 

A number of measures are used to evaluate the long-term Tb simulations 
and their associated uncertainty. Fig. 2b illustrates some of the terms used 
in this evaluation. First, we assess the quality of the deterministic Tb sim- 
ulations with the MAP parameter estimates, using the mean-square differ- 
ence (MSD [K 2 ]) between the observed and simulated long-term Tb averages 
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(Eq. 5) and standard deviations (Eq. 6) across the 24 different observations: 


MSD m = 2j- y'X m i( a MAp) - m i,o) 2 ( 5 ) 

^ 1=1 
24 

MSD s = -^2(s i (cxMAp)-Si ! o) 2 (6) 

^ i= 1 

If the modeling errors were solely due to uncertainties in the parameter val- 
ues, these metrics should be very close to zero. In practice, however, the 
metrics will substantially deviate from zero and reflect residual errors that 
cannot be explained by parameter uncertainty. The 24 differences contribut- 
ing to MSD m are illustrated as A m . in Fig. 2b. 

Secondly, we verify whether the spread in prior and posterior ensemble Tb 
simulations is in agreement with the misfit between modeled and observed 
values, in a mean-square sense. To this end, an ensemble of Tb simulations 
is generated by randomly drawing 20 samples from the prior and posterior 
parameter distributions. The misfit or skill is again defined using the mean- 
square difference (MSD [K 2 ]), but now for the ensemble means: 


MSDjn = 

kT.' 

i = 1 
24 

(mi(at) - m it 

(7) 

MSDs = 


(«*(«) - Sj,o) 

( 8 ) 


i = 1 


where “ denotes the ensemble mean. Fig. 2b illustrates the 24 differences 
contributing to MSD m as A If the uncertainties are well estimated and 
biases between observations and simulations are constrained during the cal- 
ibration, the MSD m and MSD-g metrics should match the total expected 
uncertainty ( MEnSp mi MEnSp s ), which is the sum of the Tb simulation 
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spread due to parameter uncertainty (EnSpi jmiPar , EnSpi iSiPar ) plus the resid- 
ual Tb error variance (cr 2 m , crf s ): 



(9) 


i= 1 


i= 1 



(10) 


2=1 


where cr 2 m and cr 2 s are dominated by observation, input and structural error 
after the MAP parameters values have been found. The constituent terms 
EnSp^m ,p ar and EnSpi jSjPar for each observation type i are given by: 


An illustration of EnSpi^ par is given in Fig. 2b. Again, if the uncertainties 
are well estimated, then the ratios M S D m / M EnSp m and MSDs/MEnSp s 
should be close to 1, or in other words: the “actual” ( MSD m , MSD-g) and 
“expected” ( MEnSp mi MEnSp s ) errors should be similar. These metrics 
are very similar to those used to verify the prescribed observation and simu- 
lation uncertainties in data assimilation systems (Reichle et ah, 2002) and for 
ensemble forecast verification (De Lannoy et ah, 2006). The only difference 
is that here, the mean values (i.e. the L M\ or mean, in MSD and MEnSp) 
are derived from multiple observations types (i = 1, . . . , 24), whereas in the 
earlier studies the mean was calculated in the time domain. 



(ii) 



( 12 ) 
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4. Results 


4.1. RTM- Parameters and Their Uncertainty 

In this section, we analyze the MAP values of < h >, < r > and a;, and 
their posterior uncertainty (stdv[.]). The DREAM( ZS ) scenario D a should 
be considered as benchmark in the following discussion, because of statisti- 
cal rigor of the sampled posterior (will be further discussed below). Fig. 3 
shows maps of the prior parameter values and the MAP estimates derived 
from scenario P (PSO), D and (DREAM( ZS )) (Table 1). The spatially 
averaged posterior parameter values are very similar for all 3 scenarios, with 
a microwave roughness < h > around 0.75T0.5 [-], a nadir opacity < r > 
of 0.26T0.15 [-] and a scattering albedo u of 0.09T0.07 [-], where the values 
after the ± sign measure the spatial standard deviation and reflect the vari- 
ability of the MAP parameters across the spatial domain. Note that these 
values should not be confused with uncertainty estimates. Compared to the 
prior values (Table 1 and 2), < h > has generally increased for grassland, 
< r > is smaller for forests and uj has increased for all vegetation classes 
except grassland (details per vegetation class not shown; these finding are 
similar to those of De Lannoy et al. (2013)). The spatial patterns for the 3 
scenarios are also very similar. Moreover, Fig. 3 suggests that MAP values 
derived with the PSO algorithm closely match those of DREAM( ZS ). 

Fig. 4 shows the ensemble parameter uncertainty for scenarios D and D^. 
Maps with RTM parameter uncertainty estimates for PSO (obtained as in 
De Lannoy et al. (2013)) are not shown, because they are statistically invalid 
and significantly larger than those derived with DREAM( Z s). The relative 
uncertainties for scenario D are less than 10% of the MAP parameter value 
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and substantially smaller than the spatial variability in the MAP values. For 
scenario D^, the relative uncertainties increase, with errors ranging up to 
25% of the MAP values: for < h >, the spatially averaged uncertainty is 
0.10±0.08 [-], for < r > 0.04 ±0.04 [-] and for u 0.02 ±0.02 [-], respectively. 
The uncertainty in < h > typically increases with more complex terrain and is 
smallest in the cropped region southwest of the Great Lakes. The uncertainty 
of < r > is largest in the forested Appalachian mountains where the highest 
MAP values of < r > are found. On the contrary, u is best defined in this 
area and uncertainties in uj increase in the Western dry mountain ranges. The 
< h >-values are more uncertain where either the uncertainty in u (Fig. 4e) 
or < r > (Fig. 4f) is larger. 

In summary, both DREAM( Z s) scenarios D and D a provide MAP pa- 
rameter values that are very similar and in close agreement with the PSO 
estimates. The DREAM( Z s) derived posterior parameters appear well defined 
with relative uncertainties that are less than 25% of the MAP values. It will 
be shown below that the uncertainty estimates of scenario D a are consistent 
with the sample RMSD between long-term Tb observations and simulations. 

4-2. Residual Tb Error Standard Deviation Estimation 

To analyze the effect of a m and a s in more detail, Table 2 summarizes the 
MAP parameter values and their associated uncertainties averaged over the 
entire study domain. In addition, Fig. 5 depicts the results for different veg- 
etation classes. As discussed above, scenarios D and D a return similar MAP 
RTM-parameter values, but when a m and cr s are simultaneously estimated, 
the posterior RTM-parameter uncertainty increases about 2-3 times. The 
domain-averaged values for scenario D a are a m = 3.1 K and cr s = 2.4 K 
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(Table 2), whereas scenario D uses default values of these variables of 1 K. 

The value of a m and its posterior uncertainty are largest in cropped re- 
gions (Fig. 5g) where residual Tb errors are dominated by less skillful model 
simulations. This is to be expected because irrigation is not simulated and 
the climatological LAI estimates do not account for interannual crop rota- 
tions. The parameters can not compensate for these errors, and the default 
values of a m = cr s = 1 K make scenarios D and P vulnerable to suboptimal 
solutions. For example, the relative large differences between D and D a for 
a m and a s over cropland areas increases the differences in the MAP values of 
u. For forests, a s = 1 K appears to be a good estimate (Fig. 5i) because the 
variability in Tb is expected to be low due to vegetation attenuation. Both 
the MAP values and uncertainties for cr m are always larger than those derived 
for cr s . One of the reasons for the higher a m are the opposite signs in the 
biases for the long-term averages of ascending and descending Tb, which can- 
not be mitigated with time-invariant RTM-parameters. These biases are due 
to sensor error and modeled temperature errors as discussed in De Lannoy 
et al. (2013). In a separate exercise (not shown herein), we verified that the 
a - values absorb biases in geophysical fields: by re-scaling the soil moisture 
both the RMSD (see below) and a- values are jointly reduced. 

For the simulations with prior parameters, we also calculated (i.e. not 
optimized) a m and a s as 7.5 K and 4.8 K, respectively (Table 2). Unlike 
the MAP a m - and a s - values, these prior residual a- values are dominated by 
simulation error due to suboptimal parameter values. 
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4.3. MAP Tb Performance 

Fig. 6 shows the misfit between observed and MAP simulated long-term 
Tb averages and standard deviations ( RMSD m , RMSD s: square-root of 
Eq. 5 and 6) across the 24 observations for the calibration and evaluation 
period, averaged per vegetation class. The performance skill is very similar 
for scenarios P, D and D a , which reflects that the three scenarios generate 
similar parameter estimates. The RMSD m ranges between 2 and 4.5 K dur- 
ing the calibration (Fig. 6a) and increases up to 8 K for cropland in the 
evaluation period (Fig. 6c). The RMSD s ranges between 1 and 3 K during 
calibration (Fig. 6b) and reaches values of 5 K for cropland in the evaluation 
year (Fig. 6d). Cropland has the highest errors, because of known simula- 
tion errors (see above). Note also that the RMSD m and RMSD s values of 
scenario D a during the calibration period (Fig. 6a-b) show the same pattern 
as a m and cr s in Fig. 5g and 5i. The increased errors in the evaluation period 
suggest that the calibration could benefit from climatological observations 
based on longer data records to better estimate the parameter values. 

4-4- Ensemble Tb Performance 

For DREAM( Z s), we analyze the balance between the actual Tb misfit and 
the expected uncertainty (ensemble variance) in the ensemble Tb simulations 
(20 members, as opposed to single deterministic MAP simulations above). 
The results are presented in Table 2 and Fig. 7. Table 2 shows the skill of the 
ensemble mean Tb simulations ?rq(a) and Si (a) for the calibration period in 
terms of RMSD m and RMSD-g, i.e. the square-root of Eqs. 7 and 8. These 
values are very similar to the results for the MAP simulations (section 4.3). 
For both scenarios D and D a , the RMSD m and RMSDs are respectively 
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3 K and 2.5 K, which is less than half of the actual misfit when the prior 
parameters are used. 

Table 2 also lists the square-root of the combined mean simulation and 
observation spread or expected uncertainty, i.e. RMEnSp m and RMEnSp s 
(square-root of Eqs. 9 and 10), along with the constituent terms ( RMEnSp m:Par , 
RMEnSp s , par , a m and a s ). Generally, the uncertainty associated with the 
parameter values is much smaller than the uncertainty related to other fac- 
tors, that is, RMEnSp par < cr, which is valid both when using prior and 
posterior parameter distributions. Moreover, after calibration both the o- 
and RMEnSp par -v alues are significantly reduced compared to their prior 
values. 

If the uncertainty estimates are consistent, RMSD m ~ RMEnSp m and 
RMSDj ~ RMEnSps , i.e. there should be a balance between the actual and 
expected errors (section 3.7). The domain-averaged RMSD m /RMEnSp m is 
2.7 for scenario D and 1.0 for scenario D a . Similarly, the domain- aver aged 
RMSDs/ RMEnSp s is 2.5 for scenario D and 1.0 for scenario D^. Optimal 
results are thus only found after including an estimation of a m and a s in 
scenario D a . Note that for the evaluation period (not shown), the ratios 
always exceed 1, because of an increased RMSD m and RMSDg. 

Fig. 7 shows how the ensemble spread is consistent with misfits between 
observations and simulations for scenario D a . Specifically, Fig. 7a shows the 
SMOS observed m ^ 0 and the GEOS-5 simulated rrii(ct) for ascending, El- 
polarized Tb at 6 angles for scenarios D and D^, averaged over the entire 
study domain. Fig. 7b shows the same for V-polarized Tb, and Figs. 7c 
and d provide this information for the long-term Tb standard deviations. 
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Also shown is the total ensemble simulation and observation uncertainty for 
each observation type, presented as error bars around the ensemble mean Tb 
simulations for illustration. 

The error-bars for scenario D a fully envelop the observations, whereas this 
is not the case for scenario D. Fig. 7 also explains the nature of the residual 
misfit. Except for the 57.5°-angle, the ascending ensemble mean simulations 
rrii(ct) consistently underestimate the SMOS-observed m ^ 0 for H-polarization 
and randomly deviate from the SMOS-observed mn^ 0 at V-polarization. In 
contrast, the descending simulations rrii(cx) slightly overestimate the SMOS- 
observed m^ 0 at H-polarization (see De Lannoy et al. (2013)). The SMOS- 
observed Si j0 is always larger than the simulated ^(a). This is probably 
dominated by observation noise, but could also be attributed to an under- 
estimated variability in the Tb simulations. For example, an increase in the 
RTM-parameter h not only compensates for a cold bias but simultaneously 
also reduces the Tb variability. Fig. 7 clearly illustrates why the uncertainty 
estimates obtained from scenario D a are superior. 

4-5. Convergence and Computational Cost 

The effectiveness of the posterior parameter sampling is measured by the 
convergence of the algorithms. Table 2 confirms that both the posterior 
uncertainties in the parameter estimates (stdv[.]) and the misfit between the 
simulations and observations ( RMSD ) of the long-term Tb averages and 
standard deviations are greatly reduced compared to the results with the 
prior parameter distribution. Another measure for convergence is the scale 
reduction factor, or V^R-statistic by Gelman and Rubin (1992). Values close 
to 1 are preferred, and suggest that the MCMC sampler has converged to a 
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limiting distribution. Fig. 8 shows the evolution of the convergence diagnostic 
y/R for both DREAM( Z s) scenarios. The y/R is averaged over all estimated 
parameters and across the study domain, since no obvious differences in y/R 
are found between the different vegetation classes (not shown). Initially, the 
values of y/R exhibit a lot of variation (due to random initial sample) before 
they settle down and reach values close to 1. 

Finally, we report that the derivation of the posterior distributions re- 
quires approximately 225 seconds for a single grid cell using DREAM( Z s). 
For global applications that involve 10 5 — 10 6 grid cells, posterior distribu- 
tion exploration may be too costly. Yet, if we target the MAP value only, 
PSO or DREAM( ZS ) are both viable options. 

5. Conclusions 

Accurate estimates of microwave RTM parameters for large-scale L-band 
applications are difficult to obtain. The available parameter estimates are 
generally based on small-scale field experiments and often come without any 
estimate of posterior uncertainty. This complicates radiative transfer mod- 
eling for both the forward simulation of L-band Tb over land and the re- 
trieval of soil moisture based on Tb observations. This paper expands earlier 
research reported in De Lannoy et al. (2013) to derive time-invariant RTM- 
parameters using observations of the long-term average Tb and the long-term 
Tb standard deviation obtained from SMOS data. The overall objective is to 
optimize GEOS-5 Tb simulations prior to sequential assimilation of SMOS 
or SMAP Tb data, such as planned for the SMAP L4_SM product (Reichle 
et al., 2012) and to examine the uncertainties involved in the optimization. 


25 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 


Per grid cell, 48 observations of the long-term Tb averages and standard de- 
viations were constructed for 24 different combinations of 6 incidence angles, 
2 polarizations and 2 overpass times. The differences with their respective 
long-term GEOS-5 simulations are minimized (as opposed to minimizing dif- 
ferences between Tb observations and simulations in the time domain) and 
used along with the prior parameter information to derive posterior param- 
eter estimates. 

In the present paper, the full posterior distribution of RTM-parameters 
is derived using MCMC simulation with the DREAM(zs) algorithm. To our 
knowledge, this is the first large-scale application of the DREAM( ZS ) algo- 
rithm for the estimation of RTM-parameters and their underlying uncer- 
tainty. The results serve as a benchmark to verify the results from simpler 
parameter optimization algorithms, such as for example PSO. Simple algo- 
rithms are desirable for global-scale operational applications that rely on 
evolving modeling systems in need of frequent re-calibrations. 

First, we verified that the MAP RTM-parameter values derived from 
converged posterior distributions with DREAM( Z s) can be approximated by 
a simpler optimization algorithm (PSO), which corroborates our earlier re- 
search (De Lannoy et ah, 2013). Secondly, we obtained reliable parameter 
uncertainty estimates with DREAM( Z s), which are impossible to estimate 
with PSO. The relative parameter uncertainties are generally less than 25% 
of the MAP value for < h >, < r > and u;, when including the residual 
(observation and simulation) error statistics (cr m , a s ) of the long-term Tb 
averages and standard deviations in the estimation. 

The third objective of this paper was to quantify the importance of param- 
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eter and other errors on Tb simulations. The uncertainty associated with the 
parameter values only contributes a small part to the total Tb uncertainty. 
Most of the discrepancy between Tb simulations and observations is covered 
by residual Tb errors, with MAP estimates of the standard deviations a m and 
cr s (assumed homoscedastic) around 3.1 K and 2.4 K, respectively. The prior 
estimate of 1 K was thus too low, except for cr s over forests which exhibit 
limited Tb variability due to vegetation attenuation. The largest cr m -values 
are found in cropped regions where the RMSD between Tb simulations and 
observations is also highest, due to observation errors and errors in geophys- 
ical fields (e.g. soil moisture and temperature) that constitute important 
inputs to the Tb simulations. 

The expected Tb error, i.e. the total of the MAP residual Tb error 
variance estimates (a^, <j 2 s ) and the Tb spread introduced by the posterior 
parameter uncertainties (EnSpi irriiPar , EnSpi :S:Par ), is found to be consistent 
with the actual RMSD of 3 and 2.5 K for the long-term posterior Tb aver- 
ages and standard deviations. In other words, the joint estimation of RTM- 
parameters, a m and a s with DREAM( Z s) results in a balance between actual 
and expected errors in Tb simulations, and in statistically adequate param- 
eter values and uncertainty estimates. 

In summary, the Bayesian inference of the posterior distribution of the 
RTM-parameters ensures reliable Tb simulations with GEOS-5. Further- 
more, the DREAM(zs) algorithm also reveals the importance of observation 
error and simulation error that cannot be explained by the RTM parameters. 
These error sources can be addressed using model refinement and assimilation 
of satellite-observed Tb data. 
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Appendix A. Radiative Transfer Model 


A diagnostic zero-order (tau-omega) microwave RTM is used to simulate 
L-band Tb at the top of the atmosphere (Tb toa, p [K]). The Tb toa, p at 
polarization p = [i7, V] (horizontal or vertical) is a combination of (i) soil 
emission, possibly attenuated by vegetation, (ii) vegetation emission, possibly 
reflected by the soil, and (iii) atmospheric effects: 

Tb rov,p — T s (l — r p )A p + T c (l — lu p ){1 — A p )(l + r p A p ) 

+Tb a ^ p r p A^ (AT) 

Tb toa,p Tb au,p T Gxp( T a t mp ^ T^bxov^ (A. 2) 

where Tb rov,p [K] is the top of vegetation Tb, T s [K] is the surface soil tem- 
perature, T c [K] is the canopy temperature (assumed equal to T s ), Tb a ^ p [K] 
and Tb au , p [K] are the downward and upward atmospheric radiation, A p [-] is 
the vegetation attenuation, exp(— r atrn ^ p ) [-] is the atmospheric attenuation, 
T atrn , p [-] is the atmospheric optical depth, r p [-] is the rough surface reflec- 
tivity, and uj p [-] is the scattering albedo. The atmospheric contributions 
(Tb a ^ p , Tb au , p and exp (—T atrn , p )) are described by Pellarin et al. (2003). The 
rough surface reflectivity r p [-] is derived from the smooth surface reflectivity 
R p [-] following (Choudhury et ah, 1979; Wang and Choudhury, 1981): 

r p = (Q R q + (1 - Q)Rp) exp (—h) cos Nrp (9) (A. 3) 

where Q [-] is the polarization mixing ratio and typically set to 0 for L- 
band (Kerr and Njoku, 1990), 6 [rad] is the incidence angle, h [-] is the 
roughness parameter accounting for dielectric properties that vary at the sub- 
wavelength scale, Nr p [-] is the angular dependence, and q = V for p — H 
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and vice versa. The smooth surface reflectivity R p [-] is given by the Fresnel 
equations as a function of the dielectric constant, which itself depends on soil 
moisture, temperature, texture, incidence angle and wavelength. We select 
the Wang and Schmugge (1980) soil dielectric mixing model for this study. 
The results with this model are similar to what is obtained with the Mironov 
et al. (2004) model, and both are in a better agreement with the SMOS data 
than the Dobson et al. (1985) model. We include the dependence of h on soil 
moisture (SM [m 3 .m -3 ]) through a stepwise linear expression (adapted from 
the proposed SMOS soil moisture retrieval algorithm (CESBIO et al., 2011; 
Kerr et al., 2012)): 

{ hmax if SM <= wt 

(A. 4) 

hmax + bporos^T M — wt ) if wt < SM <= pOTOS 

where poros [m 3 .m -3 ] and wt [m 3 .m -3 ] are the porosity and transition soil 
moisture, respectively. The latter is modeled as wt = 0A8.wp + 0.165 (Wang 
and Schmugge (1980)) where wp [m 3 .m -3 ] is the wilting point. 

The vegetation attenuation A p [-] is based on the Jackson and Schmugge 
(1991) vegetation opacity model: 

A p = exp( — with (A.5) 

cos 0 

Tp = b p VWC = b p LEWT LAI (A.6) 

where r p [-] is the nadir vegetation opacity, which is a function of a vegetation 
structure parameter b p [-] and the vegetation water content ( VWC ) [kg.m -2 ]. 
The latter is modeled here as the product of LAI [m 2 .m -2 ] and the leaf 
equivalent water thickness (LEWT) [kg.m -2 ]. 
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Table 1: Parameters selected for calibration in different scenarios (P, D, D a ) with an indication of the allowed parameter range 
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Table 2: Domain-averaged parameters values and their uncertainty stdv[.\ for the prior 
distributions and the posterior distributions obtained with scenarios P, D and D^. 
The bottom half of the table shows ensemble Tb prediction statistics (square-root of 
Eq. 7-8, 9-10 and 11-12), averaged across 24 long-term Tb observations and calcu- 
lated for the calibration period. Only for the prior parameters, cr m and a s are cal- 
culated assuming (a) RMEnSp m — RMSDjn and RMEnSp s = RMSDj , and (b) 
a m = and a s = sJIEEl^ - RMEnSp 2 a>par . 
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Figure 1: Study domain with indication of the dominant IGBP vegetation classes. 
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(a) 



(b) 



Figure 2: Illustration of marginal distributions for (a) RTM-parameters and (b) Tb simu- 
lations at a single grid cell. Crosses (x) indicate the MAP estimates, the vertical dashed 
lines and white box indicate the ensemble mean posterior estimate, and horizontal dotted 
arrows indicate one standard deviation uncertainty around the ensemble mean. The per- 
formance of the Tb simulations is quantified by comparing either the MAP ( rrii(a.MAP ), 
Si(cxmap)) or the ensemble mean (mj(a), s^(a)) simulations against (black dots) 24 ob- 
served values (rrii i0 , Si jQ ) with i = 1, . . . , 24. The differences A m . and A^. contribute to 
MS Dm (Eq. 5) and MSD^ (Eq. 7), respectively. 
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Figure 3: Parameter values for (left) < h >, (middle) < r >, and (right) u>, for the (top 
row) prior distribution and scenarios (second row) P, (third row) D and (fourth row) . 
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Figure 4: Uncertainty in parameter estimates for (left) < h >, (middle) < r >, and (right) 
cj, obtained with DREAM( Z g) scenario (top row) D and (bottom row) D^. 
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Figure 5: (Left) MAP parameter values and (right) uncertainties aggregated per vegetation 
class for DREAM( Z s) scenarios D and D a . Each row represents a different parameter: (a,b) 
<h>, (c,d) r, (e,f) w, (g,h) a m , (ij) a s . 
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Figure 6: RMSD in long-term Tb (a,c) average and (b,d) standard deviation during the 
(top) calibration (1 July 2011 - 1 July 2012) and (bottom) evaluation period (1 July 
2010 - 1 July 2011), using the MAP parameter values derived from PSO (scenario P) and 
DREAM( ZS ) (scenarios D and D ff ). 
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Figure 7: (a-b) Long-term average and (c-d) standard deviation, for (a-c) H- and (b-d) 
V-polarized Tb (dots) SMOS observations and (lines) ensemble simulations averaged over 
the study domain, during the calibration period (1 July 2011 - 1 July 2012) and only 
including ascending time steps. The simulations use an ensemble of parameter estimates 
derived with DREAM( Z g) scenarios (gray) D and (black) D^. The ensemble mean is shown 
by a central horizontal dash. The error bars indicate the total simulation and observation 
uncertainty and are drawn around the simulated Tb for illustration. For clarity, symbols 
are slightly offset from the nominal incidence angle. 
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Figure 8: Gelman- Rubin convergence diagnostic \/R for the two DREAM( Z g) MCMC 
simulation scenarios. The metric is averaged over all calibrated parameters, and across 
the study domain. 
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